#ifndef UHM_MAT_EXTENSION_LDL_NOPIV_HXX_
#define UHM_MAT_EXTENSION_LDL_NOPIV_HXX_

namespace uhm {
  
  // Basic LDL nopiv...
  template<class A_, class B_>
  class LDL_nopiv_ : public UHM_<A_,B_>,
                     public Direct_ {
  private:
  protected:
  public:

    LDL_nopiv_() : UHM_<A_,B_>() { }
    LDL_nopiv_(Int_ type, Int_ uplo, Int_ sym,
               Int_ mt, Int_ mb,
               Int_ rhs, Int_ blk) 
      : UHM_<A_,B_>(type, Mat_Base_::LOWER_TRIANGULAR, Mat_Base_::SYM, 
                    mt, mb, 
                    rhs, blk) { }
    LDL_nopiv_(Int_ type, 
               Int_ mt, Int_ mb, 
               Int_ rhs, 
               Int_ blk)  
      : UHM_<A_,B_>(type, Mat_Base_::LOWER_TRIANGULAR, Mat_Base_::SYM,
                    mt, mb, 
                    rhs, blk) { }
    
    LDL_nopiv_(const LDL_nopiv_ &b) : UHM_<A_,B_>(b) { }

    using UHM_<A_,B_>::set;
    using Direct_::get_leading_flop;
						

    virtual void set(Int_ type, Int_ mt, Int_ mb, Int_ rhs, Int_ blk) {
      this->set(type, Mat_Base_::LOWER_TRIANGULAR, Mat_Base_::SYM, mt, mb, rhs,blk);
    }
    
    virtual Int_ decompose(Double_ &flop) {
      Int_ mt = this->get_mt(), mb = this->get_mb();
      Bool_ complex = this->is_complex();

      flop = ( get_flop_ldl(complex, mt) +
               get_flop_trsm_lower(complex, mt, mb) +
               get_flop_syrk(complex, mt, mb) );
      return true;
    }

    virtual Int_ solve_1_x(Double_ &flop) {
      Int_ mt = this->get_mt(), mb = this->get_mb(), rhs = this->get_rhs();
      Bool_ complex = this->is_complex();

      flop = ( get_flop_trsm_lower(complex, mt,rhs) +
               get_flop_gemm(complex, mb, rhs, mt) );
      return true;
    }

    virtual Int_ solve_2_x(Double_ &flop) {
      Int_ mt = this->get_mt(), mb = this->get_mb(), rhs = this->get_rhs();
      Bool_ complex = this->is_complex();

      flop = ( get_flop_trsm_upper(complex, rhs,mt) +
               get_flop_gemm(complex, mt, rhs, mb) );
      return true;
    }

    virtual Int_ update() {

      if (this->get_mt() && this->get_mb()) {
        // ** Make ATL = LU ... may be this is better
        this->_ATL.hermitianize();
        apply_diag_scal( Mat_Base_::LEFT, Mat_Base_::NO_CONJUGATE,
                         Mat_Base_::UPPER_TRIANGULAR, 
                         this->_ATL );
        trsm( Mat_Base_::RIGHT, Mat_Base_::UPPER_TRIANGULAR, 
              Mat_Base_::NO_TRANSPOSE, Mat_Base_::NONUNIT_DIAG,
              1, this->_ATL, this->_ABL );

        /*
          trsm( Mat_Base_::RIGHT, Mat_Base_::LOWER_TRIANGULAR, 
          Mat_Base_::TRANSPOSE, Mat_Base_::UNIT_DIAG,
          1, this->_ATL, this->_ABL );
          apply_diag_inv_scal( Mat_Base_::RIGHT, Mat_Base_::NO_CONJUGATE,
          this->_ATL, this->_ABL );
        */

        // ATR = DTL*ABL^T 
        this->_ATR.create_flat_buffer();
        copyt( this->_ABL, this->_ATR );
        apply_diag_scal( Mat_Base_::LEFT, Mat_Base_::NO_CONJUGATE, 
                         this->_ATL, this->_ATR );

#pragma omp taskwait

        gemm_l( Mat_Base_::NO_TRANSPOSE, Mat_Base_::NO_TRANSPOSE,
                -1, this->_ABL, this->_ATR, 1, this->_ABR );

        this->_ATR.free_flat_buffer();
      }
      return true;
    }
    
    virtual Int_ decompose() {

      if (this->get_mt()) 
        ldl_nopiv( Mat_Base_::LOWER_TRIANGULAR, this->_ATL );

      this->update();

      return true;
    }

    virtual Int_ solve_1_x() {
      if (this->get_mt()) {
        trsm( Mat_Base_::LEFT, Mat_Base_::LOWER_TRIANGULAR,
              Mat_Base_::NO_TRANSPOSE, Mat_Base_::UNIT_DIAG,
              1, this->_ATL, this->_XT );
      }
      if (this->get_mt() && this->get_mb()) {
        gemm( Mat_Base_::NO_TRANSPOSE, Mat_Base_::NO_TRANSPOSE,
              -1, this->_ABL, this->_XT, 1, this->_XB );
      }
      return true;
    }

    virtual Int_ solve_2_x() {
      if (this->get_mb() && this->get_mt()) {

        this->_ATR.create_flat_buffer();

        copyt( this->_ABL, this->_ATR );
        apply_diag_scal( Mat_Base_::LEFT, Mat_Base_::NO_CONJUGATE,
                         this->_ATL, this->_ATR );

        gemm( Mat_Base_::NO_TRANSPOSE, Mat_Base_::NO_TRANSPOSE,
              -1, this->_ATR, this->_XB, 1, this->_XT );

        this->_ATR.free_flat_buffer();
      }

      if (this->get_mt()) {
        trsm( Mat_Base_::LEFT, Mat_Base_::UPPER_TRIANGULAR,
              Mat_Base_::NO_TRANSPOSE, Mat_Base_::NONUNIT_DIAG, 
              1, this->_ATL,  this->_XT );
      }
      return true;
    }

    virtual Int_ solve_1_r() {
      if (this->get_mt()) {
        trsm( Mat_Base_::LEFT, Mat_Base_::LOWER_TRIANGULAR,
              Mat_Base_::NO_TRANSPOSE, Mat_Base_::UNIT_DIAG,
              1, this->_ATL, this->_RT );
      }
      if (this->get_mt() && this->get_mb()) {
        gemm( Mat_Base_::NO_TRANSPOSE, Mat_Base_::NO_TRANSPOSE,
              -1, this->_ABL, this->_RT, 1, this->_RB );
      }
      return true;
    }

    virtual Int_ solve_2_r() {
      if (this->get_mb() && this->get_mt()) {

        this->_ATR.create_flat_buffer();

        copyt( this->_ABL, this->_ATR );
        apply_diag_scal( Mat_Base_::LEFT, Mat_Base_::NO_CONJUGATE,
                         this->_ATL, this->_ATR );

        gemm( Mat_Base_::NO_TRANSPOSE, Mat_Base_::NO_TRANSPOSE,
              -1, this->_ATR, this->_RB, 1, this->_RT );

        this->_ATR.free_flat_buffer();
      }

      if (this->get_mt()) {
        trsm( Mat_Base_::LEFT, Mat_Base_::UPPER_TRIANGULAR,
              Mat_Base_::NO_TRANSPOSE, Mat_Base_::NONUNIT_DIAG, 
              1, this->_ATL,  this->_RT );
      }
      return true;
    }

  };

}

#endif
